TIFR preprint no. TIFR/TH/ 00-34 



accepted for publication in J. Stat. Phys. 



Continuously varying exponents in a sandpile model 
with dissipation near surface 

S. Liibeck 1 and D. Dhar 2 

1 Theoretische Tieftemperaturphysik, Gerhard-Mercator-Universitat Duisburg, 

Lotharstr. 1, 47048 Duisburg, Germany 

2 Department of Theoretical Physics, Tata Institute of Fundamental Research, 

Homi Bhabha Road, MUMBAI, 400005, India 

February 1, 2008 
Abstract 

We consider the directed Abelian sandpile model in the presence of sink 
sites whose density f t at depth t below the top surface varies as ct~ x . For 
X > 1 the disorder is irrelevant. For \ < 1, it is relevant and the model is no 
longer critical for any nonzero c. For % = 1 the exponents of the avalanche 
distributions depend continuously on the amplitude c of the disorder. We 
calculate this dependence exactly, and verify the results with simulations. 
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1 Introduction 

An important question in the general area of self-organized criticality Q is the 
universality of critical exponents. Do small perturbations of the local evolution 
rules leave the critical exponents unchanged? How do perturbations of the system 
at the boundary affect the critical steady state? These are the questions we study 
in this paper, and construct a simple analytically tractable model where changes 
near the boundary lead to continuously varying critical exponents. 

The influence of boundary perturbation in nonequilibrium systems is more com- 
plicated than in equilibrium systems ||. In the latter, away from phase transition 
points, the bulk is largely unaffected by the boundary except that in case of spon- 
taneous symmetry breaking, the boundary conditions will remove the degeneracy 
between different phases, and some finite-size effects. Near critical phase transi- 
tions, the critical behavior near the boundary is substantially modified, and con- 
stitutes a different class of (surface) critical phenomena, characterized by a new 
set of (surface) exponents ||. In the so-called 'ordinary transitions', the surface 
phase-transition is driven by the bulk transition in the sense that the transition 
occurs at the critical temperature corresponding to the bulk, and its position or 



the critical exponents are not affected by local changes in the Hamiltonian near 
the surface. 

In self-organized critical systems, the correlation length is, by definition, always 
large, and the influence of perturbations at the boundary is felt deep inside. In 
the well-known Bak-Tang-Wiesenfeld (BTW) sandpile model of self-organized crit- 
icality even the existence of a critical steady state depends on the presence of 
a surface where particles can leave the system. Thus, in such systems, boundary 
perturbations affect the overall behavior of the system to a much greater extent. As 
the steady state is critical, various measured average quantities in the steady state 
near a surface will differ from the bulk values by an amount which will decrease as 
a power law of the distance from the surface. 

For example, in the BTW model in two dimensions, Brankov et al. || have 
shown that the mean height of the pile at a site at distance t from the surface is 
less than the mean value in the bulk by an amount proportional to t~ 2 . In a one 
dimensional sandpile model with threshold dissipation, Ali found that the density 
is larger near the boundary || and the excess density varies as 1/t. We would 
like to understand the effect of such extended perturbations on other properties 
of the system. For example, one may ask if such a density profile significantly 
affects the probability of extinction of avalanches that reach the boundary. One 
could even ask if it determines the avalanche exponents In this context, it 
seems worthwhile to understand the properties of self-organized critical systems 
with 'artificially produced' power-law surface perturbations. 

In this paper, we consider the directed Abelian sandpile model in the presence 
of some sink sites, whose density varies with the distance from the surface with a 
power x- We are able to determine the x-dependence of the avalanche exponents 
exactly. Our analysis shows that the disorder is irrelevant if x > 1> an d is relevant 
for x < 1. For x = 1> the perturbation is marginal, and the critical exponents 
depend continuously on the disorder strength. In the case of equilibrium critical 
phenomena, for Ising model with deviation of the coupling strength from the critical 
value varying with distance from surface, a similar behavior is seen H, and may 
be understood in terms of a simple renormalization argument ||. 

In fact precisely this question comes up in the sandpile model studied by Tadic 
and Dhar (TD) |10|| . This model (its definition is given later in the paper) belongs 
to the universality class of directed models with stochastic toppling rule, with no 
multiple topplings For this model, TD presented some intuitive arguments to 
determine the exact critical exponents for a general dimension d in terms of the 
exponents of the directed percolation problem. The arguments, involving deter- 
mining the statistics of clusters in a directed percolation model with space- varying 
density, are plausible, but not rigorous. In the model studied here, our exact an- 
swer coincides with that obtained based on the same reasoning, thus provides some 
additional evidence in the latter's favor. 

In the next section we define our model. We determine the disorder-averaged 
two-point correlation function of toppling events exactly in section ||. This shows 
the qualitative difference in the behaviors of the avalanche distribution for x > 1 
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Figure 1: The lattice structure of the directed Abelian sandpile model. Sand-grains 
are added at the top and removed from the bottom. Periodic boundary conditions 
are used in the horizontal direction. The full circles correspond to toppled lattice 
sites and the open circles mark sinks where sand-grains leave the system in the 
bulk. 



and x < 1- We give heuristic arguments to determine the avalanche exponents 
as a function of the amplitude c in the marginal case % = 1. These results are 
confirmed by numerical investigations. A somewhat more careful derivation of 
the avalanche exponents requires the consideration of the three-point correlation 
function of toppling events which is presented in section |j. The connection to the 
argument used in flfC]], and some generalizations to other systems are discussed in 
the concluding section. 



2 Definition of the model 

We consider a two dimensional lattice in the shape of an open cylinder (see Fig. [l]). 
Each site is labelled by two integer coordinates (x, t), where x labels the coordinate 
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in the horizontal direction, and t denotes the depth of the layer (0 < t < T). It is 
convenient to work with only the set of sites for which x + 1 is even (one of the two 
sublattices into which this lattice may be decomposed). We denote the number of 
sites in each horizontal row by L. To each lattice site (x,t), a non-negative integer 
variable h(x, t) is associated, called the height of the sandpile at that site. A sand- 
grain is added at a randomly chosen site x on the top layer (t = 0). If the height 
of a site (x,t) exceeds the critical value h c = 2, i.e, if h(x,t) > 2, it is unstable 
and it topples. On toppling at a site, its height is decreased by 2, and one grain is 
transferred to each of the two downward neighbors (x =L l,t + 1). In this way the 
transferred sand-grains may activate the two downward sites and an avalanche of 
relaxation events may take place. Some typical avalanches are shown in Fig. 

The size of an avalanche may be measured in terms by its downward length 
(called its duration t), or the number of toppled sites (called its mass m). In 
the critical steady state, the probability distribution of avalanche duration and 
avalanche mass has power- law tails [jl2| , i.e., 

3 

Piob{duration > t} ~ t~ Tt+1 with r t = — , (1) 

Prob{mass > m} ~ m~ Tm+1 with r m = — . (2) 

In this paper, we consider the effect of introducing surface disorder in the model. 
We do this by declaring a small fraction of sites as sink sites. Any particle falling 
onto a sink site is lost from the system. The sink sites are uncorrelated with each 
other. We assume that the concentration of sink sites ft at the depth t decreases 
as a power of t: 



ft = c t~ x , for t » 1. (3) 

The case where the sink concentration is uniform was studied earlier by Tadic et 
al. |l3j . They showed numerically that any finite concentration of sinks destroys the 
criticality of the steady state. Theiler |l4|] showed that in the mean-field theory 
a uniform sink density / induces an exponential cutoff in the depth reached by 
avalanches which varies as f^ 1 for small /. 

In the next sections we investigate the two-point and three-point correlation 
functions which allow to derive the avalanche exponents exactly. 



3 The two-point correlation function 

The characterization of the steady state is very simple in our model. From the 
general theory of Abelian sandpile models IB], it is easy to see that all stable 
configurations are recurrent, and occur with equal probability in the steady state. 
Thus, in the steady state, each non-sink site is empty or occupied with equal 
probability, independent of the occupation of other sites. 
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Let us now consider the propagation of avalanches in the system. Let G(x,t\x' ,t') 
be the probability in the steady state that a toppling will occur at (x, t) if a particle 
is added at (x',t'). If the site (x,t) is a sink site, we define G(x,t\x' ,t') to be zero. 
In the steady state, at non-sink sites, the average number of particles coming into 
(x, t) must equal the number of leaving particles. Thus, the two-point correlation 
function G(x,t\x't') satisfies the equation 

G(x,t + l\x',t') = -rj(x,t + l) G(x - l,t\x',t') + G(x + l,t\x',t') , for t>t', 

(4) 

where rj(x, t + 1) is zero if the site is a sink site, and one otherwise. Averaging 
over the configuration of sink sites gives the disorder-averaged two-point correlation 
function Ct(x, t+1 \ x', t'). Using (rj(x,t+l)) = (1 — /t+i) and the fact that r](x, 
is independent of G{x ± l,t\x',t'), it is easy to see that G(x,t + l\x',t') satisfies 
exactly the difference equation 



G(x,t+l\x',t') = i(l-/ m ) [G(x-l,t\x',t') + G(x + l,t\x',t') 



for t > t'. 

The disorder-averaged correlation function G is translationally invariant in the 
x-direction, i.e., 

G(x + x ,t\x' + x ,t') = G(x,t\x',t'), (6) 

for all ^o- We use Eq.((5|) as a recursion equation and determine G(x,t\0, 0) for all 
t, given its initial value at t = 

G(x,0|0,0) = i(l-/ )^). (7) 

It is easy to show that the disorder-average function G is given by 

G(x,t\0,Q) =V t G o (x,t\0,0) (8) 
where Go is the propagator for the problem with no sinks, i.e., 

G (MO,0) = 2-'- 1 (9) 

and where Vt is defined to be the product 

Vt = II C 1 - /;)■ ( 10 ) 

3=0 

For the case with no disorder, it is known that Prob{duration > t} and 
Go(0, t|0,0) both vary with the same power of t | |12| |. If we make the plausible 
assumption that this continues to hold in the presence of sinks, we expect that the 
probability distribution of the avalanche duration is of the form 

Prob{duration > t} ~ t^l 2 V t (11) 
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Figure 2: Snapshot of five avalanches for various values of x an d c. In all cases the 
avalanche duration is t = 500. In order to display the details the avalanches are 
stretched by the factor 2 in the horizontal direction. 



since the propagator Go is known to decrease asymptotically as t 1//2 for large 

* El- 

Let us consider some limiting cases. In the case of homogeneous disorder ft = f 
we get that Vt ~ (1 — /)*, and the probability that an avalanche reaches depth t 
decreases exponentially with t. This is in agreement with the numerical results of 
& 

In the case that the disorder density ft decreases as a power-law according to 
Eq. @ one has to distinguish three cases: 

i) For x > 1) the disorder is irrelevant. Vt tends to a finite constant value for 
t — > oo, and the asymptotic behavior of the disorder averaged correlation function 
G is the same as that of Go, i.e., the disorder does not affect the scaling behavior. 
A typical avalanche for x = 2 is shown in Fig. While branchings can occur, most 
of the longer avalanches have only one prominent branch. 

ii) For x < 1) the disorder is relevant, and leads to loss of criticality in the 
system. Vt tends to zero as exp (— where K is a constant which depends 
on x an d c. Thus the function G decreases to zero with the distance as a stretched 
exponential. In contrast to the case x > 1 the avalanche propagation for x < 1 is 
characterized by many branchings, though only one branch survives for long (see 
Fig- D- 

Hi) For x = 1) the disorder is marginal and the avalanche distribution exponent 
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Figure 3: The probability distribution Prob{duration > t} for \ = 1 an d various 
values of c for a system of size L x L. The dashed lines correspond to the derived 
scaling behavior [Eq. (]lT|)]. The probability distributions Pr ob {duration > t} are 
averaged over 50 different representations of the disorder (hole-configuration) and 
2 10 non-zero-avalanches for each disorder representation. 

varies continuously with the parameter c. Here, Vt tends to zero for large t as t~ c . 
From Eq. (||), we see that G decreases as i _1 / 2_c and thus we find 

r t = \ + c. (12) 

The avalanche distribution Prob{duration > t} still exhibits a power-law decay, 
but the exponent depends continuously on the amplitude c. This behavior has 
been verified in our simulations. Figure || displays the probability distribution 
Prob{duration > t} for three different values of c. For sufficiently large durations 
the asymptotic behavior agrees with Eq. (|i~2|). 

Let us now consider the distribution of masses of avalanches. We see from the 
explicit form of the avalanche propagator G that while its dependence on t is mod- 
ified by the disorder, its dependence on the transverse coordinate x is unaffected. 
Hence, even in the presence of disorder, the average transverse size of an avalanche 
at depth t scales as t 1 / 2 . As each site in an avalanche topples only once, the average 
width of an avalanche cluster varies as t 1//2 and the total number of toppled sites 
m scales as i 3 / 2 . Then Prob{mass > m} ~ m~ Tm+l with 

4 2 

Tm = 3 + 3 c - ( 13 ) 
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Figure 4: Finite-size scaling analysis of the mass probability distribution 
Probjmass > m} for \ = 1> three different values of c and for various system 
sizes L £ {256, 512, 8192}. To avoid an overlap the curves for c = 0.5 and 
c = 0.25 are shifted in the lower-left direction. 



In our simulations, we used lattices of size L x T with T = L. In this case, 
the transverse size is effectively infinite as even the largest avalanches are not 
affected by the finite transverse size. Thus the finite-size effects are governed by 
the longitudinal size (T) alone, and the mass distribution is expected to fulfill the 
finite-size scaling ansatz fllq] 

Prob(mass > m|T) ~ T' 13 f(T~ u m). (14) 

If finite-size scaling works a data collapse of the distributions for different system 
sizes LxT (with L = T) is obtained for v = 3/2 and j3 = 1/2+c. The corresponding 
scaling plots for three different values of c are shown in Fig. ||. In all cases we 
observed good data collapse which confirms the above analysis. 



4 The three-point correlation function 

In the previous section we obtained the avalanche exponents by using some scaling 
arguments, which are not careful in distinguishing the probability that a given site 
at depth t is toppled in an avalanche from the probability that at least one site at 
depth t is toppled. We now present a cleaner derivation of the above results. 
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Let us denote the number of sites toppling in the layer t by N(t). Based on 
the experience with the model without disorder, we make the assumption that the 
distribution of avalanches in this model can be described by the scaling function 

P(N(t)=n) ~ At- a F(n/t b ) (15) 

for n^O and where the scaling function J-{x) decreases fast for large x, and A is 
a constant. From this scaling form, it follows that (N(t)) ~ t 2b ~ a , and (N 2 (t)) ~ 
i 3b-a . Thus knowing the t-dependence of (N(t)) and (N 2 (t)), we can determine 
the exponents a and b. Clearly, using Eq.(@) we get 

(N(t)) = G(x,t|0,0)= lp t . (16) 

We would now like to calculate (N 2 (t)). Let G (xi,x 2 ,t) be the probability that 
both sites x\ and x 2 at depth t topple when a particle is added at x = on the top 
layer t = 0, averaged over all configurations of the disorder. Then, clearly, we have 

(^ 2 (*)) = J2T,G i3 \xi,x 2 ,t). (17) 

XI X2 

(3) 

For t > 0, the three-point correlation function G satisfies the equation 

G i3 \x 1 ,x 2 ,t) = (l-f t ) 2 J2 \G i3) (x 1 + e l ,x 2 + e 2 ,t-l) (18) 

ei,e 2 =±l 

for x\ ^ x 2 . We have to supplement this equation with the boundary condition for 
x\ = x 2 : 

G i3 \ Xl ,x u t) = G( Xl ,t). (19) 
As for the case with no disorder, Eq.(|i~8|) is solved by the ansatz 

t 

G (3) (x!, x 2 , t) = E E /(2/> u ) G(xi,t\y, u) G(x 2 ,t\y, u) (20) 
n=0 y 



which satisfies Eq.(18) automatically for x\ ^ x 2 . We choose the function f(y,u) 
so that the equation holds also if x% = x 2 . Putting x\ = x 2 = x in the previous 
equation, we get 

t 

G(x,t) = ^2Y,f(y^ u )G 2 (x,t\y,u) (21) 
u=o y 

This is a linear equation connecting f(y,u) to its values at sites higher up. Hence 
these can be solved recursively. The calculation of (N 2 (t)) can be further simplified. 
We define the function 

=£/(!/,«), (22) 
y 
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and note that 



Y,G(x,t\y,u) =±Vt/V u -i, 

X 

where we define V—i = 1. Then using Eqs.(17j20|), we get 

(N 2 {t)) = -V 2 J2 V-\F(u). 



(23) 



(24) 



u=0 



The recursion equation for F(u) is also easily deduced from the recursion equation 
for the function f(y, u). We note that J2 X G ( x i AVi u ) ls independent of y and only 
depends on t and u. Let us write 



J2G 2 (x,t\y,u) = K(t,u). 



(25) 



Then summing Eq.(21) over x, we get 

1 



■Pt = Y,F( U )K(t, 

u=0 



(26) 



The effect of the disorder is just to introduce an overall rescaling factor to K(t, u) 
in the absence of disorder: 



K(t,u) =V 2 t Vz\ £ G 2 Q (x,t\y,u). (27) 

X 

This implies that the asymptotic behavior of F is given by 

F(u) ~ u- c - 1/2 . (28) 

Substituting in Eq.(|24|), we see that 

(N 2 (t)} ~ t 1 ' 2 ' . (29) 

The obtained scaling behavior of (N(t)) and (N 2 (t)) implies that the exponents 
of the assumed scaling form [Eq. (|l5|)l are given by a = 1 + c, and b = 1/2. This 
justifies the earlier heuristic calculation of the exponents. 

It is straight-forward to extend this treatment to higher dimensions. We only 
mention the result here: in two (i.e. 2 + 1 dimensions) and higher dimensions the 
avalanche exponents are given by 

r t = 2 + c, and r m = (3 + c)/2. (30) 



As in the case with no disorder, there are logarithmic correction terms to the 
exponent values given above in the two-dimensional case. 
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5 Discussion 



Let us comment briefly on the relationship of this model to that studied earlier by 
TD. Their stochastic variant of the directed sandpile model is defined on the same 
lattice as presented in Fig. [l], and particles are added on the top layer at random, 
and removed from below. The difference from the model studied here is that there 
is no quenched randomness, and the stochastic toppling rules are defined as follows: 
If a particle is added to a site, only then is the site examined for toppling. If the 
number of particles exceeds 2, with a non-1 probability p a toppling occurs and 
one particle is transferred to each of the two downward sites, and with probability 
(1 — p) nothing happens. While there is particle conservation in the TD model, so 
far as the evolution of a single avalanche is considered, the sites where no toppling 
occurs when particles are added, act as sink sites. The density of sink sites is 
higher near the top in both models. The randomness in the TD model provided 
by the stochastic toppling rules is 'annealed', unlike the quenched randomness in 
the model studied here. 



In [10], it was shown that this model has a steady state for p > p*, where p* 
is the critical probability for directed site percolation on a square lattice. In the 
steady state, there is a small fraction of sites having number of particles 2, 3, 4, . . ., 
but the average density is finite, and is a function of the depth of the layer. 

TD argued that a density less than the bulk density locally corresponds to a 
finite correlation length which is spatially varying. They argued if the medium has 
an inverse correlation length n(t), which is small, and varies slowly as a function 
of t, then the survival probability of a disturbance up to depth t is given by 

Prob(t) ~ Probo(t) exp (- J* n(t')dt^j . (31) 

The first factor is what we should get if n(t) is zero everywhere. The second 
term is a simple generalization of the exp(— nt) factor when varies with t. 
We note that the expression is non-perturbative, in that the small parameter k is 
exponentiated. More importantly, k is usually a nontrivial power of the starting 
perturbation 5p(t) in the density (in the TD case, re(i) ~ [6p(t)] v , where v is the 
longitudinal correlation length of directed percolation problem). 

In the TD model, n(t') varies as c/t' , and hence the argument of the exponential 
varies as — c logt, for large t. This immediately gives the c-dependence of the 
power-law behavior of Prob(f). 

Using explicit calculation of correlation functions for small finite t, it is easy 
to verify that Eq. ( |3l|) is not exactly valid for the TD case. However, for the 
subsequent derivation the exponents of the TD model to be valid, the ansatz has 
to be at least asymptotically exact. As a proper treatment of the correlations in 
the steady state of TD model seems very difficult, one may try to find other simpler 
models where the validity of the ansatz can be demonstrated. 

The steady state in our model is very simple, and thus the model is a good 
testing-ground. It is gratifying that, in the present model, the ansatz does become 
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exact, and the two-point correlation function satisfies Eq. (|8|) exactly. Thus we 
have shown that this approximation is exact for our model. 

The approximation is in the spirit of the familiar semiclassical Wentzel-Kra- 



mers-Brillouin- Jeffreys approximation in quantum mechanics [17|. We note that 
the mathematical mechanism giving rise to continuous variation of critical ex- 
ponents with amplitude of perturbation applies to a much wider class of prob- 
lems. Consider, for example, the generalization of the equilibrium problem treated 
in |H involving extended surface defects to general p-vector spins {S(n)} in d- 
dimensional half-space. Lattice sites are labelled by a d-dimensional integer vector 
n = (ni, ri2, • • • , rid) with n\ > 0. The Hamiltonian of the system is 

U = J(n,n')S(n) • S(n') (32) 

where the summation is over all nearest neighbor pairs of sites n, n'. The coupling 
constant J(n, n') is assumed to depend only on ri\ (we may assume that n\ < n^) 
as 

J(n,n') = J*-CV (33) 

where J* is the value of coupling in the bulk, and C > 0. 

As mentioned earlier, the simple renormalization argument due to Burkhardt 
H shows that for \ > V^j where v is the (bulk) correlation length exponent, this 
perturbation is irrelevant, and for \ < 1/u, it is relevant. In the marginal case 
X = we get the inverse correlation length ft(ni) ~ C v jn\. Let us consider the 
spin-spin correlation function (S(n) • S(n')). For simplicity, we may take n and 
n' to differ only in the first coordinate. Then, it is easy to write an expression 
for this correlation function using the ansatz Eq. (j3ll). The mathematical problem 
then is the same as discussed in this paper, and one finds that in the marginal case 
X = and the critical exponent corresponding to the decay of this correlation 
function varies linearly with C v . 

Of course, the dramatic effect of the boundary-perturbations on the exponents 
of the avalanche size distribution is due to the fact that all avalanches are initiated 
at the top boundary in our model. If we add particles at a fixed distance to from the 
boundary, it is easy to see that large avalanches of duration much greater than to 
have the same asymptotic behavior as the case to = 0. However, shorter avalanches 
find significantly fewer defects than if they had been added at the top. If we added 
particles everywhere in the bulk, the avalanche size distribution is averaged over 
different values of to, and in the thermodynamical limit is the same as the case 
without defects. 

We thank M. Barma for a critical reading of the manuscript. S. Liibeck would 
like to thank V. B. Priezzhev and A. Hucht for useful discussions. 
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